其他
Biostar课程：5、6
#为作者的注释
#*为译者的注释
Lecture 5 - 编写可重复使用的脚本。获取子序列。
# 通过locus号来获取一条序列。
efetch -db nucleotide -id KM233090 -format fasta > KM233090.fa# 这条序列也可以通过accession号来获取。
efetch -db nucleotide -id 667853062 -format fasta > 667853062.fa# 可以一次性获取多个id下的序列。
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta > multi.fa# 查看文件有多少行。
wc -l multi.fa# 数有多少个头文件。
cat multi.fa | grep ">" | wc -l# 查看头文件。
cat multi.fa | grep ">"# 我们也可以从数据中只获取一个子集。
efetch -db nucleotide -id KM233090,KM233066,KM233113.1 -format fasta -seq_start 1 -seq_stop 10 > multi.fa# 我们开始来获取所有的埃博拉病毒基因组的开头。
#*这条命令获取的是每条序列的第1-10个碱基。
esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10 > starts.fa# 每次都手动输入命令太费事，是时候使点黑技能来解决这个问题了。开始写shell脚本吧。
# shell脚本就是很多命令的集合，是可以重复运行的。
# 把下边的这些行放到一个文件中，通常以.sh为后缀。我将其命名为get_seq.sh
# esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start 1 -seq_stop 10
#现在，你可以把一个脚本里的所有命令运行了。
bash get_seq.sh > starts.fa# 能从外部指定这个脚本的某个部分，会更好。
# 两个$符号将被来自命令的参数取代。 详情参见下边的内含文件。
bash get_seq.sh 1 10 > starts.fa# 现在，只把这个文件中的序列提取出来。
# -A 1 是指把匹配行及其下一行打印出来。
# -v 是指把不匹配的行打印出来。
# 不要直接一下子输入全部语句，而要每输入一个语句，检查结果，再输入第二个，如此反复。
cat starts.fa | grep ">" -A 1 | grep -v ">" | grep -v "-"我们第一版的脚本：
#!/bin/sh## 从所有的序列中把核苷酸子序列抓取出来。# 存放在PRJNA257197 这个项目（project）里。## 预期有start和stop这两个参数。# 这个设置使脚本停止在错误和未定义变量上。set -ue# 通过Entrez Direct访问NCBI。esearch -db nucleotide -query PRJNA257197 | efetch -db nucleotide -format fasta -seq_start $1 -seq_stop $2我们前面脚本的升级版。当你加上函数时，生成不同的命令就变得更为容易了：
# 从所有的序列中把核苷酸子序列抓取出来。# 存放在PRJNA257197 这个项目（project）里。# 预期有start和stop这两个参数。# 这个设置使脚本停止在错误和未定义变量上。set -ue
function esearch_ebola() { # 搜索ebola（埃博拉）的核苷酸序列。
esearch -db nucleotide -query PRJNA257197
}
function efetch_limits() { # 你可以把变量表示什么变得更明确。
# 对于更为复杂的命令，这应该会帮助。
START=$1
END=$2
# 运行efetch.
efetch -db nucleotide -format fasta -seq_start $START -seq_stop $END}# 把两个命令链起来获取数据。# 请留意着两个参数是如何传递到函数中的。esearch_ebola | efetch_limits $1 $2Lecture 6 - 编码和fastq格式
# 使用python后跟着-c这个命令来运行字符串。
python -c 'print "Hello World"'# 字符的序数（整数）。
python -c 'print ord("A")'# 整数的字符(字母代码) 。
python -c "print chr(65)"# 找到字符对应的Sanger（+33）编码质量值
python -c 'print ord("I") - 33'# 找到字符对应的Illumina 1.5 (+64)编码质量值
python -c 'print ord("I") - 64'# 找到40这个值对应的Sanger (+33)编码质量字符
python -c 'print chr(40 + 33)'# 找到40这个值对应的Illumina 1.5 (+64)编码质量字符
python -c 'print chr(40 + 64)'#* 关于这个质量值和字符的对应关系，在《小白生信学习记3》有详细介绍。
# 从这个课程的网址中获取一个测试数据集。
curl -O http://www.personal.psu.edu/iua1/courses/files/2014/illumina-data.fq.gz#* 原文是curl -O
#* 这个地址是当时上课的学生用的。我们是不能这么用的。
# 你可以解压这个gzip文件 illumina-data.fq.gz ，也可以不解压，并用gzcat (Mac OSX)或者zcat (Linux)来打开。当你有很多的数据集时，牺牲读取速度来节省存储也是值得的。
time gzcat illumina-data.fq.gz > tmp# 解压数据，同时保留原压缩文件。
# 没有 -c 的话，你可以无需重定向一个输出，它会在当前文件夹下解压。
gunzip -c illumina-data.fq.gz > illumina-data.fq
time cat illumina-data.fq > tmp# 查看FASTQ记录 (Mac OSX)
gzcat illumina-data.fq.gz | head -4# 查看FASTQ记录 (Linux)
zcat illumina-data.fq.gz | head -4# 拆包一个SRA数据和机器下来的原始数据集
# 此处会生成2个文件。
fastq-dump --split-files SRR1553605#*这个数据之前是下载过的。详情参照上两节的翻译。
# 显示其中一个fastq的记录。
cat SRR1553605_1.fastq | head -4# 简单粗暴地检查一下质量，在前10,000行里是否有连续的###字符。
gzcat illumina-data.fq.gz | head -10000 | grep '###' | wc -lcat SRR1553605_1.fastq | head -10000 | grep '###' | wc -l# 哪个看起来要好些？
# 只看前面/最后的几行的质量是有些冒险的，因为数据的排序可能不是随机的。通常我们会生成一个随机样本（我们在后边会看到如何去实现）。
# 我们可以在数据里搜索pattern并着色。
# 这些都是起始密码子吗？是否存在pattern TATA?
gzcat illumina-data.fq.gz | head -10000 | grep 'ATG' --color=always | head
gzcat illumina-data.fq.gz | head -10000 | grep 'TATA' --color=always | head# 有看起来像Illumina接头: GATCGGA的序列吗？
gzcat illumina-data.fq.gz | head -10000 | grep 'GATCGGA' --color=always | head翻译：生物女博士
注：本系列课程翻译已获得原作者授权。
欢迎关注我们